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Abstract 

I present a hybrid-like two-step algorithm, which combines a microcanonical update of a 
spin system using demons, with a multicanonical demon refresh. The algorithm is free from 
the supercritical slowing down, suffered by canonical methods: the exponential increase of the 
tunnelling time between the metastable states in the first-order phase transitions, when the 
volume of the system is increased. The demons act as a buffer between the multicanonical heat 
bath and the spin system, allowing the spin system to be updated with any microcanonical 
demon procedure, including cluster methods. The cluster algorithm is demonstrated with the 
2-dimensional 7-state Potts model, using volumes up to 128^. The tunnelling time is found 
to increase as L^'^^, where L is the linear dimension of the system. 
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1 Introduction 



If a system has a first-order phase transition, at the transition temperature it can exist in a 
mixed state, where two different bulk phases are separated by an interface. The free energy 
carried by the interface is the interface tension a. Because at the transition temperature the 
free energy densities of the bulk phases are equal, the free energy of the mixed state is higher 
than either of the pure phases by an amount of Fs = crA, where A is the area of the interface. 
In numerical simulations using the canonical ensemble at the transition temperature, the 
configurations containing the mixed phase are suppressed by the Boltzmann factor e"^*/^. 
When the volume of the system is increased, the suppression of the mixed state increases 
exponentially with the area of the interface - hence also the time it takes for the system to 
tunnel from one pure phase to another increases exponentially. 

Recently, Berg and Neuhaus introduced a powerful new method, the multicanonical 
algorithm, which avoids the exponential slowing down by artificially enhancing the probability 
of the configurations with an interface. The tunnelling time increases only polynomially with 
the increasing linear size L. In this method the individual spin updates generally depend 
on the total energy of the system in a non-linear fashion. This effectively prevents the use 
of vector or parallel coding (except if one runs many lattices in parallel) and cluster update 
algorithms. 

In the microcanonical demon Monte Carlo approach, developed a decade ago by M. Creutz 
||2|, one uses additional variables, demons, to transfer energy around the system. The total 
energy of the system plus the demons is absolutely conserved. However, by periodically 
updating the demon energies according to the Boltzmann distribution the method reduces to 
the canonical procedure. 

In this work I present an algorithm which combines these two methods: first, a system is 
updated microcanonically with a set of demons, and second, the demons are refreshed with 
a multicanonical heat bath. The demons act as a buffer, isolating the actual system from 
the multicanonical heat bath, thus enabling one to choose an optimal microcanonical update 
step for a particular problem. The update can be a highly vectorizable and parallelizable 
local update or a cluster update, or some combination of these. As an example, I apply the 
cluster version of the algorithm to the 2-dimensional 7-state Potts model. 

The Potts models have become standard tools for high-precision Monte Carlo studies 
of the first-order phase transitions. The 2-dimensional g-state (2dg's) Potts model ^ is 
defined by the partition function 




(1) 
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E{s) = Y,{l-6{s„s,)), Si = l,...,q. (2) 

When g > 4, the transition is of first order. The infinite volume transition temperature is 
Pc = l/7c = log(l + ^/q)- Beside the fact that many infinite volume quantities are exactly 
known, rigorous finite-size scaling (FSS) predictions by C. Borgs et al. offer a quantitative 
method for studying the approach to the asymptotic regime. I chose the 2d7s Potts model for 
comparison with the recent standard multicanonical calculation by W. Janke et al. and 
the canonical one by A. Billoire et al. [^. The lattice sizes in this work were V = = 20^, 
32^, 64? and 128^; for the largest volume, two separate simulations were performed. 

This article is divided into three parts: in the first section I discuss how the standard 
multicanonical approach is generalized to the two-step demon algorithm, and how it can be 
used to obtain an estimate of the density of states, and, through this, of all thermodynamical 
quantities. The next section describes the actual update algorithms: the microcanonical 
cluster update and the multicanonical demon refresh. I present a method which enables 
the 'slow' part of the multicanonical update to be performed in oc VV steps. The results 
of the 2d7s Potts model simulations are reported in the third section. The tunnelling time 
was found to increase like L^-^'^^^\ which is better than the standard multicanonical method 
|-^2.65^_ Where appropriate, the thermodynamical measurements of the 128^ lattices fully 
agree with the rigorous FSS predictions of ref. [^, and all the measurements are consistent 
with the multicanonical MC data of ref. ||^]. 



2 The Multicanonical Demon Algorithm 

Close to the transition temperature the canonical probability distribution pp{E) develops a 
double-peak structure, and one definition for the transition temperature = itself 
is the temperature when the two peaks have equal height (fig. 1). Following refs. Q, 
the peak locations are denoted by and E2 ^ and the probability density is normalized to 
Pp{E^) = pj^{E2) = 1. Denoting the minimum of pj^ between the peaks by p^i^, the interface 
tension is [^]: 

a = -lhn^^^. (3) 

Because of the periodic boundary conditions, the configurations corresponding to the mini- 
mum of the probability distribution have two interfaces separated by ~ L/2; hence the factor 
2 in eq. (^). In the following the L-dependence of the above quantities is mostly suppressed 
from the notation. 

In the standard 'direct' multicanonical method, one usually aims at a roughly constant 
probability density in the domain between the peaks: pwiE) = 1, Ei < E < E2- This can 
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be achieved by substituting the Boltzmann weight with a weight function W{E): 

pw(.E)^ns(,E)e-'^^^\ (4) 

where ns{E) is the number of states at energy E. The requirement that the probabihty 
density is constant imphes that W{E) oc logn5'(-E) = S{E) when Ei < E < E2; and 
W{E) = (3cE + const, otherwise. Because ns{E) is what we are trying to compute in the first 
place, we have to use an approximate W{E) instead - obtained, for example, with finite-size 
scaling, canonical simulations, or previous multicanonical simulations. The measured pw is 
then reweighted with e^^^^ to produce ns{E), from which all quantities of interest can be 
calculated - also the improved W{E). 

Now we want to apply multicanonical ideas to a system consisting of the original spin 
system and a system of demons. Heuristically, it is clear that the spin system in eq. (Q) can 
be substituted with any other system, also with this composite system. Denoting the weight 
function with G, the probability density can be written as a joint distribution in the spin 
system energy Es and the demon energy E^: 

Pg{Es, Ed) oc ns{Es)nD{ED) e'^^^^) , (5) 

where ns{Es) and noiEj:,) are the spin system and the demon density of states, respectively, 
and Et = Es + -Ed is the total energy. In the most general case the weight G is a function 
of both Es and Et, but when Et is fixed, we want pc to reduce to the microcanonical 
distribution ns{Es)nD{ET — Es), implying that G can only be a function of Et- When Es 
is fixed, pg reduces to the multicanonical distribution for the demons: nTt^Eo) e~^^^^~^^°\ 
Thus the probability distribution @ is preserved in a generic two-step process of a micro- 
canonical spin update and a multicanonical demon update, provided that both of the update 
steps separately satisfy detailed balance. Note that if G{E) = f3E, both the demon and spin 
system will have canonical distributions. 

After the Monte Carlo simulation has produced a sample of the distribution (|5|), ns can 
be solved from it. Although not strictly necessary, it is advantageous to use the known demon 
density of states: with Nti demons with discrete energy states 0, 1, . . ., 

,p , {ND-1 + En)l 
{Nn-l)lEnl " 

Because eq. (|5|) is valid separately for each Ed, ns can be expressed as a linear combination 

^ /I Pg{Es,Ed) . , . 

ns{Es) oc Ae^;Es iEn)e-G^^T) ' 2^ ^^^^^s = 1 ■ (7) 
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The multipliers A^j^-Es should be chosen to minimize the error in n^. For each Eg separately, 
the uncertainty in the (measured) pc is given by Spc = \Jpg/Nes) where pc is the probability 
distribution in the limit of infinite statistics and N^g is the number of measurements with this 
Es- (This is not true for the full distribution pg{Es, Ed) because of the correlations between 
successive measurements; however, for fixed Es, each measurement of Ed is completely 
independent, as explained below). Minimizing the resulting error in eq. we finally obtain 

"'^^'^ Eed nliED)e-OiEr) ' Ps{Es)^Y:pg{Es,Ed). (8) 

Note that the final result depends only on the spin system distribution ps{Es), not on the 
demon energy distribution. The distribution ps corresponds to the standard multicanonical 
distribution eq. @ with the weight 

^-wiEs) ^ ^ ^^^Ed) e-^(^«+^^) . (9) 

Ed 

Quite generally, if we want to simulate a system with a non-canonical weight function W{E), 
then by inverting eq. (^) we obtain the corresponding G{E). With the two-step update, the 
function G will produce exactly the same distribution as W with a direct update. 

Instead of aiming at a fiat distribution of Es, it is now more natural to try to 'fiatten' the 
Et distribution. Then, the optimal G{Et) equals to St{Et), the entropy of the total system. 
In the actual runs described here the initial estimate of G{Et) was obtained (for the lattices 
< 64^) from short runs with G{Et) = PEt; for the 128^ lattice, finit e-size scaling was used 
to scale up the function used in the 64^ simulation. The 64^ and 128^ lattices required one 
further refinement run to obtain the final weight function; in the end two different weights 
were used for the 128^ lattice. To simplify the calculation, Et was restricted in the range 
< Et < EJ^^^, where the limits were chosen such that the expectation values of Es, as 
a function of Et, bracket the peak locations Ei and E2 (fig. 1): 

{Es){Ef-) <E,<E2< {Es){En- (10) 

The functional form of G{Et) is not crucial, as long as it is accurate enough; here a continuous 
piecewise linear form was used. The required accuracy increases with the volume of the system 
- the weight G{Et) is an extensive quantity (or rather G{Et) — PcEt oc L), but if G{Et) is 
'wrong' by an amount of, say, log 2, the probability pg{Et) will be changed by a factor of 2. 

The left part of fig. 2 shows the joint distribution from the 64^ lattice, using Es and Et as 
independent variables. As can be seen, Es and Et follow each other very closely - the length 
and the width of the ridge behave like and L, respectively. The demon energy varies only 
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very little, meaning that the microcanonical temperature dS{E) / dE is almost constant, as it 
should be in the phase coexistence region. The right part is the same distribution, 'canonized' 
by reweighting it with qG{Et)-PcEt ^ where f3c is the transition temperature for this lattice 
(table I). 

An interesting variation of the algorithm can be obtained by restricting Et to a discrete 
set of values £^0) • • • > E^^, sufficiently dense so that the neighboring Es distributions have 
large enough overlaps. This is a microcanonical version of the 'simulated tempering' method, 
presented by Marinari and Parisi^]. 



3 The Update Algorithms 

3.1 The Microcanonical Cluster Update 

In the simulations reported here one update cycle consists of one microcanonical spin + demon 
update sweep followed by an energy measurement and a multicanonical demon update. As 
the microcanonical step I used the Swendsen-Wang variation of the microcanonical cluster 
algorithm presented recently by M. Creutz [10|. As opposed to the standard procedure, the 



demons are located on the links of the lattice, instead of being connected to the spins. A 
link is activated only if the spins have the same value at each end of it and the demon does 
not carry enough energy to frustrate the link. Clusters are grown over the activated links 
and each cluster is flipped to a random spin value. Finally, the demon energy is increased 
or decreased by the amount of the corresponding link energy decrease or increase. On a 
d-dimensional lattice, this method requires d y~V demons. 

After each update cycle the demon locations are shuffled. If the shuffling were not per- 
formed, one would construct exactly the same clusters during the next update cycle - as- 
suming that the demon refresh, described in the next chapter, is also skipped. The shuffling 
does not need to be perfect: here it was done with random offsets and step lengths when 
picking the demons from the demon array. The actual cluster search was performed with the 



Hoshen-Kopelman cluster-finding algorithm [11 1 



Note that it is also possible to perform local updates with the same demons; the demons 
can be left on the links or moved to the spins - in the latter case only half (or l/d) of the 
demons are used during one sweep. The probability distribution (^) is unaffected, as long as 
the total number of demons remains the same. 



3.2 The Multicanonical Demon Update 

The multicanonical demon refresh is greatly facilitated by the fact that each demon is an 
independent degree of freedom and the demon density of states is known. The most straight- 
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L 


iterations 


TL 


a 


20 


2 500 000 


320(5) 


0.0189(3) 


32 


5 000 000 


821(15) 


0.0169(2) 


64 


6 000 000 


2700(81) 


0.0147(4) 


128a 


9 000 000 


10720(520) 


0.01302(17) 


1286 


6 000 000 


10520(620) 


0.01306(21) 



Table 1: The tunnelling time and the interface tension from the 2d7s Potts model simulations. 



forward way to perform the demon update is to touch each demon with the multicanonical 
heat bath: the demon i is assigned a new value with the weight exp 

For a continuous demon energy this would be the best method, even though this is a non- 
vectorizable process with ~ Nd steps. However, since now the demon energy is discrete, 
there exists a method which enables a major part of the demon update to be performed in 
~ V^D (non-vectorizable) steps. First, a new total demon energy S^f™ is calculated with a 
global heat-bath update: let j; be a random number from even distribution between and 1; 
now the new demon energy is the smallest satisfying 

-Y.E"=onD{E")e-0{E"+Es)^ ^ ^ 

where G{Et) = oo, when Et < i?™™ or Ej- > Ej^^^. This guarantees that the demon energy 
at fixed Es is free from autocorrelations, justifying the use of the multinomial distribution 
prior to eq. (P). A new demon state with energy E'ff^ can then be constructed from the old 
one by adding or subtracting energy from randomly selected demons. However, care has to 
be taken to ensure proper counting of states: energy is added or subtracted unit by unit, and 
the demon to be changed is chosen according to the respective probabilities 



To prove that eq. ( [1^ ) is correct, it is sufficient to show that starting from the state E^ = 0, 
produces with equal probability all the states with the same energy. Let us construct a 
specific demon state oj with an energy E^^■. by repeatedly applying the probability of a 
particular sequence of additions {i} leading to this state becomes 

E^ = Ni + 2N2 + W^ + ... , (14) 
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where A^e is the number of demons with energy e. Because the numbers N^, are characteristic 
of the state uj and not of the sequence {i}, all the sequences leading to this state have the 
same probability. Then the probability p^^ of producing the state uj is obtained by multiplying 
Pjj} with the number of sequences N^^^^ (= number of permutations) leading to this state, 

giving = l/n£){E^) [eq. @]. This is just what we want - all the states with the same 
energy are equally probable, and the sum of the probabilities is 1. It is obvious that p_ is 
the probability of the inverse process. Because the old demon state can be understood as an 
intermediate step in the energy addition/subtraction sequence, we can start building the new 
state directly from it. 

Remember that during one update cycle only additions or only subtractions are performed. 
Because the average energy change in one cycle is ~ VV, only as many demons need to be 
updated; the demon shuffling after each cluster update takes care of proper mixing. By 
careful use of auxiliary arrays, also the integrations in eq. ([l^) can be performed in \fV 
steps, but this was not fully implemented in the simulations. 

An important technical question is how to implement the demon selection according to 
the probabilities of eq. (|l2|). Two simple methods were tested: first, one can employ an 
appropriate accept/reject step at a random demon location - this is a true oc V-^D method, 
but the poor acceptance rate (8% for the 2d7s Potts model) makes this rather slow in practice. 
In the second method one constructs a pointer array a, which has an entry for the location 
of each unit of the demon energy. The array o has Ed elements, of which E'^j^^ are pointing to 
the demon i. When Ed Ed + 1, the demon is chosen as follows: one generates a random 
integer i in the range \\-^Nd + Ed\i \ii < Nd-, the demon i is selected, and if i > Nd, one 
chooses the demon Oj. By construction, this gives just the correct probability p\. The demon 
address is then added to the array a, and Ed is increased. The process is repeated until i?]~f™ 
is reached. The subtraction of the energy is performed analogously; only now does one select 
a random pointer Oj, decrease the energy of the demon Oj, and set = 0. The bottleneck in 
this method is the initial generation of the pointer array, which requires ~ Nd steps; however, 
it is a vectorizable process, and since Ed has to be measured anyway, it can be performed 
with little overhead. In the tests the pointer array method was found to be 5-8 times faster 
than the accept/reject procedure, so it was chosen for the actual simulations. On the 128^ 
lattice, the whole measurement and demon update cycle used less than 6% of the CPU time, 
the rest being taken up by the cluster operations, which contain the only non-vectorizable 
oc V loop. Using a Cray X-MP, one full cycle took ^ 3.9 /xs per spin. 
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4 Results 



The performance of the algorithm is best measured by the tunnehing time, which is defined 
as in ref. 0: four times the tunnelhng time, Atl, is the average number of updates for the 
system to get from to E2 and back. This definition gives comparable autocorrelation time 
definition. The times are listed in table |l|, and shown in fig. 3. The fit to the three largest 
lattices gives tl = 1.49(17) x L^-^'^(^) . This scales better than the standard multicanonical 
method |], which has tl = 0.082(17) x l2-65(5), fact, this is even better than the optimal 
scaling given by the random-walk picture: in the first-order transitions, the energy gap 
behaves as V, but the width of the system energy distribution with fixed total energy increases 
only like Assuming an ideal update algorithm, this is also the average change in the 

system energy during one sweep. The system has to random- walk across the gap in order 
to tunnel from one phase to another, and this takes (gap/step)^ ~ y = sweeps. The 
discrepancy is largely due to the shift of Ei /V and E2 /V (see figs. 1 and 8), which have a 
finite size dependence ~ 1/L. If we ignore these and calculate tl from tunnellings between 
the infinite volume energy density values ei = 0.4453 . . ., 62 = 0.7986 . . ., the tunnelling times 
of the smaller lattices become shorter and we obtain a scaling law tl = 0.66(7) x L^-^'^(^\ 
which is compatible with the random-walk limit. However, I chose to present the former 
result in order to enable the comparison with the previous calculation. The quality of the 
scaling fit is only = 6.2/2 d.o.f., and since the 'physical' correlation length is of order ~ 30 
iP, it is quite plausible that the true scaling law will be different. The scaling functions are 
plotted in fig. 3. 

In order to compare with the canonical update algorithm, I utilized the results of the 
2d7s Potts model simulations by A. Billoire et al. Their work has high-statistics data 
from 5 different volumes between 16^ and 64^. I made a finite-size fit to the autocorrelation 
times with the heuristic function tl = aL" e^'^^ , where the interface tension a had the 
fixed value 0.01174 (see the next paragraph). The parameter values given by the fit are 
a = 1.01(15) and a = 2.31(4), with = 4.2/3 d.o.f.. (Without the exponential factor, the 
best power-law fit gives = 56/3 d.o.f.) The function is plotted in fig. 3. The simulations in 
ref. 1^ were performed with the one-hit Metropolis algorithm, which has a notoriously long 
autocorrelation time; this was balanced by a fast multispin coding. With a heat-bath update 
the autocorrelation time could be ~ 5-8 times shorter, bringing the low end of the line close 
to the level of the multicanonical lines. On the 128^ lattice the multicanonical cluster method 
is ~ 3 times faster than the standard multicanonical, which again is ~ 5-40 times faster than 
the canonical method. 

The interface tension was measured with the Binder method Q , eq. . The minima and 
maxima of the canonical probability distributions were found by fitting a parabola close to 
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L 


/3(equal height) 


/?(C'max) 


/3(-Bmm) 


/?( equal weight) 


20 


1.28474(13) 


1.28443(12) 


1.28444(13) 


1.2939(3) 


32 


1.28976(7) 


1.28953(7) 


1.28776(7) 


1.29379(7) 


64 


1.29241(4) 


1.29235(3) 


1.29194(3) 


1.29360(4) 


128a 


1.293251(15) 


1.293234(15) 


1.293140(15) 


1.293567(16) 


1285 


1.293242(19) 


1.293227(19) 


1.293133(19) 


1.293560(19) 



Table 2: Measured pseudotransition temperatures. 



the extrema; the results depend on the fitting range only very weakly, when the range is large 
enough. All the error analysis was done by jackknifing the data to 50 sets. The measured 
values of a are shown in fig. 4, together with the measurements of ref. Q. Using a common 
FSS fit of the form ^ 

- ^ logPmin = ^ + X 

to the three largest volumes, we obtain the result a = 0.01174(19) with c = 0.169(11). The 
result agrees well with ref. (c = 0.0121(5)), but is seven standard deviations off from the 
exact infinite- volume value a = 0.010396..., recently calculated by Borgs and Janke [|l2[. 
The cited errors are only statistical, and the large difference between the values is obviously 
caused by the violations of the FSS formula (16). This is supported by the missing 'flat 
bottom' around the Pmin' ^^i^ part corresponds to the variations of the distance between 
the two interfaces, and the absence of the flat part implies that the interaction between the 
two interfaces is still non-negligible. 

The value of a is at least a factor of 6 smaller than the results of refs. []T3| , p!4| , obtained 
with an unrelated method. In these calculations the interface was stabilized by using different 
temperatures on the different sides of the lattice, thus enforcing one side into the disordered 
phase and the other into the ordered one. The interface tension was measured as a function 
of the temperature difference, and the final answer was obtained by extrapolating the results 
to the transition temperature. This method has also been applied to A'^^ = 2 pure gauge QCD 
p^ , 17 1, and the results agree with the values obtained with the Binder method ||^, |T^. The 



apparent failure of this method in the case of the 2d7s Potts model is probably due to the 
too strong pinning effect caused by the temperature difference. In 2d, the amplitude of the 
interface fluctuations is large, ~ ^jL/aT. These fluctuations are strongly suppressed by the 
temperature difference, which should then be very small; however, in that case the two-phase 
configuration would be lost, unless extremely large volumes are used. 

Finally, let us compare transition temperature measurements to the exact finite-size ex- 
pansions [^, 20 1 . Common definitions for the pseudotransition temperature are the locations 
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of the maximum of the heat capacity C = jV {{E"^) — (E)'^) and the minimum of the 
Binder parameter Vl = ^{1 — {E'^) / {E'^)'^), where E' = E — 2V in order to comply with the 
definition of energy used in [^,0]. The differences of the known infinite- volume transition 
temperature and the measured values are plotted in fig. 5; also shown are the exact lowest 
order FSS corrections. On the 128^ lattices the FSS correction is within the error bars of the 
measured values, whereas the 64? lattice is still off by ~ 2a. Figure 6 shows the behaviour 
of Vl as a function of (3. 

Still another way to find the transition temperature is the 'equal weight' method [19|, 



where is defined as the temperature when the relative probabilistic weights of the ordered 
and disordered states are q = 7 and 1, respectively: 

q = Wo/Wd = E Pp^i^)/ E P/3^(^)' (17) 
E<E' E>E' 

where E' is the energy at the minimum of pp at the temperature when the two peaks have 
equal height. The measurements of f3c are shown in fig. 7. The FSS corrections in this case 
are only exponential. An FSS ansatz of the form = + a e~^^ was fitted to the data, 
with the results a = 0.0012(12), b = 0.05(3), and (3^ = 1.293562(14). The fit gives exactly 
the correct (3c', however, is almost completely determined by the 128^ data, and the effect 
of the values of a and b to the value of is negligible. The various transition temperature 
measurements are listed in table ^. 

The locations of the maxima of the canonical probability distribution p/^ are shown in 
fig. 8. The FSS fits to the three largest volumes give the infinite-volume result ci oo — 
0.4421(24) and 62,00 = 0.7981(17), which agree fairly well with the exact values 0.44539 . . . 
and 0.79867. . .. The latent heat is the difference of these two values: Ae = 0.3559(29) (exact 
0.35327...). 

Even though the difference between the heat capacity of the disordered phase (Cd) and of 
the ordered phase (Co) is exactly known, the actual pure phase values are not. An estimate 
can be obtained by employing the FSS relation of refs. [^, 0]: Co = Cx^max — F(As/2)^ + 
0.0038 . . . + 0{V~^), where As = (3c^E/V is the entropy difference between the two phases. 
The result is shown in fig. 9; the FSS fit yields an estimate of Co = 44.4 it 2.2. Because of the 
subtraction of two terms of order V, the errors grow rapidly when the volume is increased. 
Again, this result is consistent with that of ref. 47.5 it 2.5. 

5 Conclusions 

I have presented a new hybrid-like algorithm, which combines a microcanonical spin system 
update with demons, and a multicanonical demon update. Like the direct multicanonical 
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method, this algorithm does not suffer from the exponential slowing down in the first-order 
phase transitions. In the 2d7s Potts model simulations, the tunnelling time was found to in- 
crease like tl ~ L^'^^(^^ with lattices up to = 128^. Where appropriate, the measurements 
were compared with the analytical finite-size scaling formulas by Borgs et al. Q; within the 
statistical errors, the 128^ lattice was found to be in complete agreement with the order 1/V 
FSS predictions. Also, all results are fully compatible with those of ref. 0. However, the 
common FSS ansatz for the interface tension, = o"oo + c/L, fails to produce the correct 
infinite-volume value. This is probably due to the interactions between two interfaces, which 
are ignored by the FSS ansatz. Clearly, a better FSS function is needed. 

The main advantage of the multicanonical demon algorithm is that it offers various possi- 
bilities in choosing the algorithm. In addition to simply using either local or cluster updates, 
one can also adjust the number of the demons and the number of the microcanonical updates 
before each multicanonical step. For example, if one has a very fast local microcanonical 
algorithm, it might be preferable to interleave many microcanonical updates for each demon 
refresh, and to use a large number of demons {Nd > V) in order to allow large fluctuations in 
the system energy during the microcanonical phase - the demon refresh can still be performed 
in Nu X fast -|- \/Nd x slow operations. The method can also be generalized to magnetic 
transitions by using demons carrying magnetization; in this case the cluster algorithm cannot 
be used. 
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Figure captions 

Fig. 1. The canonical probability distributions pp^{E) obtained from the multicanonical 
distributions with eq. (^. 

Fig. 2. The multicanonical joint probability distribution pg{Es, Et), measured from the 64^ 
7s Potts model simulation (left) and the canonical distribution corresponding to /? = 1.29241 
(right). The canonical distribution is obtained from the multicanonical one by reweighting it 
with e^(^^)-/^^^. 

Fig. 3. The tunnelling time as a function of the lattice size L on a log- log scale with scaling 
laws. Solid line (a): this work, oc L^'®-^; dash-dotted line (b): multicanonical simulation 
by Janke et al. 0, L^'^^; dashed line (c): canonical simulation by Billoire et al. L^'^ e^"^ . 

Fig. 4. The interface tension as a function of the lattice size. The black circles correspond 
to this work, and white squares to measurements by Janke et al. 0. The line is a fit to 3 
largest volumes, and the white circle is the infinite volume extrapolation. The arrow is the 
exact infinite volume result |12[ |. 

Fig. 5. The difference between the infinite volume transition temperature and the mea- 
sured transition temperature (5^ , plotted on a log-log scale. The black circles are the measured 
locations of the maxima of the heat capacity (Cmax) and the white circles the locations of 
the minima of the Binder parameter (Vz^^niin)- The solid line and the dashed line are the 
exact order 1/V expansions. For clarity, the data of the second 128^ lattice has been shifted 
slightly to the right. 

Fig. 6. The Binder parameter V^. The dashed line shows the exactly known infinite volume 
limit. 

Fig. 7. The transition temperature determined by the equal weight method, and the FSS fit 
to the data. The dashed line is the exact infinite- volume /3c- 

Fig. 8. The locations of the maxima of the canonical probability distribution p^, together 
with the linear fit to three largest volumes. 

Fig. 9. The maximum of the heat capacity with the leading finite size correction subtracted. 



13 



